In [1]:
import plotly.express as px
import plotly.graph_objects as go
import plotly.io as pio #graph output
from plotly.subplots import make_subplots
from ipywidgets import interact, widgets
import pandas as pd
pd.options.plotting.backend = "plotly"
import numpy as np
In [2]:
def tribo_cond_df(file_no):
    if file_no<=9:
        filename = 'Results/CSV/PUEG4/PV_Limite_00{}.csv'.format(file_no)
    else:
        filename = 'Results/CSV/PUEG4/PV_Limite_0{}.csv'.format(file_no)
    df1 = pd.read_csv(filename, header = 17, skiprows = [21])
    if file_no<=9:
        filename = 'Results/Chronoamperometry/PUEG4_recip_0{}.csv'.format(file_no)
    else:
        filename = 'Results/Chronoamperometry/PUEG4_recip_{}.csv'.format(file_no)
    df2 = pd.read_csv(filename, header = 4)

    df1 = df1[['T','COF']]

    df2.columns = ['T','i']

    df1['COF1000'] = df1.COF.rolling(1000).mean()
    df2['i10'] = df2.i.rolling(10).mean()
    return df1,df2,file_no
def plot_tribo_cond(df1,df2,file_no):
    fig = make_subplots(specs=[[{"secondary_y": True}]])
    fig.add_trace(go.Scattergl(
        x=df2['T'], y=df2['i'],
        name='Corrente',
        mode='lines',
        marker_color='rgba(152, 0, 0, .8)'
    ),
                  secondary_y=False,
                 )
    fig.add_trace(go.Scattergl(
        x=df1['T'], y=df1['COF'],
        name='COF',
        mode='lines',
        marker_color='rgba(0, 0, 152, .8)'
    ),
                 secondary_y=True)
    fig.update_xaxes(range = [705,706])
    fig.update_yaxes(title_text = 'Corrente (µA)', color = 'rgba(152, 0, 0, .8)', range=[-0.0006,0.0006], secondary_y = False)
    fig.update_yaxes(title_text = 'Coeficiente de Atrito (-)', color = 'rgba(0, 0, 152, .8)', range=[-0.5,0.5], secondary_y = True)
    fig.update_layout(title_text = 'PUEG4 - Ensaio {}'.format(file_no))
    #df = pd.merge(df1[['T','COF']],df2[['T','i']], on = 'T')
#     fig2 = px.density_heatmap(x=df.loc[df['T']>10,"COF"],
#                              y=df.loc[df['T']>10,"i"],
#                              marginal_x="histogram",
#                              marginal_y="histogram"
#                             )
    #fig3 = px.scatter(df2[df2['T']>10], x = 'T', y='i', trendline = 'lowess')
    return fig
In [3]:
df1, df2, file_no = tribo_cond_df(2)
In [4]:
fig1,fig2,fig3 = plot_tribo_cond(df1,df2,file_no)
In [5]:
df2['i20'] = df2.i.rolling(20).mean()
df2['i20s'] = df2.i.rolling(20).std()
In [6]:
df2['Cond'] = df2['i10']>1.96*df2['i20s']
In [7]:
df2[df2['T']>5].plot(kind = 'scatter', x = 'T', y = 'i', color = 'Cond')
In [8]:
df1,df2,file_no=tribo_cond_df(2)
fig1 = plot_tribo_cond(df1,df2,file_no)
In [9]:
fig1
In [10]:
fig1.write_html('tribocond.html')
In [11]:
df = pd.DataFrame()
forca = {1:2.50,
         2:2.50,
         3:5.00,
         4:2.50,
         5:5.17,
         6:5.17,
         7:8.00,
         8:8.00,
         9:8.00,
         10:5.17,
         11:5.00,
         12:8.00,
         13:2.50,
         14:5.00,
         15:5.00,
         16:5.17,
        }
freq = {1:2,
         2:4,
         3:4,
         4:4,
         5:3,
         6:3,
         7:2,
         8:4,
         9:2,
         10:3,
         11:4,
         12:4,
         13:2,
         14:2,
         15:2,
         16:3,
        }

for k in range(1,17):
    df1, df2, file_no = tribo_cond_df(k)
    df2['i20'] = df2.i.rolling(20).mean()
    df2['i20s'] = df2.i.rolling(20).std()
    df2['Cond'] = df2['i20']>1.96*df2['i20s']
    df = df.append({'Ensaio':k,
                    'Tempo Conduzindo %': df2.loc[df2['Cond'],'T'].count()/df2['T'].count(),
                    'COF Médio': df1.loc[df1['T']>0.75*df1['T'].max(),'COF'].mean(),
                    'Corrente Média':df2.loc[df2['T']>30,'i'].mean(),
                    'Frequência': freq[k],
                    'Força Normal': forca[k]                    
                      }, ignore_index=True)
In [12]:
df['Resistência Média']=3000000/df['Corrente Média']
In [13]:
df
Out[13]:
COF Médio Corrente Média Ensaio Força Normal Frequência Tempo Conduzindo % Resistência Média
0 0.174521 0.000056 1.0 2.50 2.0 0.301436 5.389928e+10
1 0.172148 0.000069 2.0 2.50 4.0 0.144069 4.329921e+10
2 0.142043 0.000013 3.0 5.00 4.0 0.043483 2.244868e+11
3 0.169514 0.000034 4.0 2.50 4.0 0.113368 8.918589e+10
4 0.156568 0.000030 5.0 5.17 3.0 0.026562 1.007814e+11
5 0.168654 0.000045 6.0 5.17 3.0 0.016862 6.605898e+10
6 0.159651 0.000052 7.0 8.00 2.0 0.246972 5.735881e+10
7 0.137771 0.000040 8.0 8.00 4.0 0.139546 7.583892e+10
8 0.145512 0.000037 9.0 8.00 2.0 0.142342 8.025238e+10
9 0.154457 0.000047 10.0 5.17 3.0 0.024623 6.418894e+10
10 0.157316 0.000040 11.0 5.00 4.0 0.060208 7.481072e+10
11 0.146632 0.000048 12.0 8.00 4.0 0.043625 6.303415e+10
12 0.159981 0.000042 13.0 2.50 2.0 0.289088 7.059272e+10
13 0.152554 0.000041 14.0 5.00 2.0 0.153916 7.339449e+10
14 0.144294 0.000036 15.0 5.00 2.0 0.202605 8.383175e+10
15 0.145532 0.000030 16.0 5.17 3.0 0.052670 1.002929e+11
In [14]:
s = []
gl = 0
for f in df['Força Normal'].unique():
    for fr in df['Frequência'].unique():
        mask = (df['Força Normal']==f)&(df['Frequência']==fr)
        if df.loc[mask,'Corrente Média'].size>0:
            s.append(df.loc[mask,'Corrente Média'].std())
            gl += df.loc[mask,'Corrente Média'].size-1
s = sum(s)/gl
In [15]:
fig = px.scatter_matrix(df[df['Resistência Média']<2e11],
                  dimensions=["COF Médio", "Resistência Média", "Tempo Conduzindo %", "Força Normal", 'Frequência'],
                       color = 'Corrente Média')
fig.update_traces(diagonal_visible=False)
fig.update_layout(
    title='Condutividade in-loco',
    dragmode='select',
    width=990,
    height=990,
    hovermode='closest',
)
fig.show()
In [16]:
fig = px.scatter_3d(df,
                    x='Força Normal',
                    y='Frequência',
                    size ='Tempo Conduzindo %',
                    z = 'COF Médio',
                    color = 'Resistência Média')
In [17]:
fig.update_layout(coloraxis_colorbar_tickformat = 's',
                  scene_camera = dict(
                      up=dict(x=0, y=0, z=1),
                      center=dict(x=0, y=0, z=0),
                      eye=dict(x=1.35, y=1.35, z=1.35)
                  )
                 )
In [18]:
fig = px.scatter_3d(df, x='Força Normal',
                    y='Frequência',
                    z='Resistência Média',
                    color = 'COF Médio',
                    size='Tempo Conduzindo %',
                    size_max=30,
                    opacity = 0.7,
                   )
fig.update_layout(scene_zaxis_tickformat = 's',
                 scene_camera = dict(
                      up=dict(x=0, y=0, z=1),
                      center=dict(x=0, y=0, z=0),
                      eye=dict(x=1.35, y=1.35, z=1.35)
                  ),
                 )
In [19]:
df['Intercept']=1
df['Ff']= df['Força Normal']*df['Frequência']
df['F²']= df['Força Normal']**2
df['f²']= df['Frequência']**2
P = df.columns.to_list()
In [20]:
P
Out[20]:
['COF Médio',
 'Corrente Média',
 'Ensaio',
 'Força Normal',
 'Frequência',
 'Tempo Conduzindo %',
 'Resistência Média',
 'Intercept',
 'Ff',
 'F²',
 'f²']
In [21]:
parametros = [P[7], #intercept
              #P[3], # força
              #P[4], # freq
              #P[8],
              #P[9],
              P[0]
              #P[10]
             ]
resposta = P[1]
In [22]:
X = df[parametros].to_numpy() #matriz de parâmetros X
Y = df[resposta].to_numpy().T #matriz de respostas Y

b = np.dot(np.dot(np.linalg.inv(np.dot(X.T,X)),X.T),Y) #cálculo do vetor de parâmetros b 
Vb = (np.linalg.inv(np.dot(X.T,X)))*s**2 #variância/covariância dos parâmetros b, com base na variância calculada na primeira parte do programa
Vb = np.diag(Vb).T #organizada como o vetor b (exclui covariancias)
In [23]:
b
Out[23]:
array([-6.06094952e-05,  6.55217674e-04])
In [24]:
Vb
Out[24]:
array([1.04548804e-09, 4.30488130e-08])
In [25]:
def anova(data,entradas,resp,coef,coef_var,matriz,nome):
    from scipy.stats import f
    #cálculo das variáveis para a ANOVA
    #respostas calculadas com o modelo
    data = data.sort_values(by = ['Força Normal','Frequência','COF Médio'])
    y_hat = 0
    for i in range(len(coef)):
        param = data[entradas]
        y_hat = y_hat + (coef.item(i)*param[param.columns[i]])

    #respostas médias para as repetiçoes de ensaios
    y_bar = []
    m = 0
    #for eg in np.sort(data['EG'].unique()):
    for forca in np.sort(data['Força Normal'].unique()):
        for freq in np.sort(data['Frequência'].unique()):
            for cof in np.sort(data['COF Médio'].unique()):
                mask = (data['Força Normal']==forca)&(data['Frequência']==freq)&(data['COF Médio']==cof)
                if data.loc[mask,resp].size>0:
                    m = m+1
                    for i in range(data.loc[mask,resp].size):
                        y_bar.append(data.loc[mask,resp].mean())


    p = matriz.shape[1] #número de parâmetros do modelo
    n = data[resp].count() #número total de observações

    SQR = ((y_hat-data[resp].mean())**2).sum() #soma quadrática da regressão
    GLR = p-1 #graus de liberdade da regressão
    MQR = SQR/GLR #média quadrática da regressão

    SQr = ((data[resp]-y_hat)**2).sum() #soma quadrática dos resíduos
    GLr = n-p #graus de liberdade dos resíduos
    MQr = SQr/GLr #média quadrática dos resíduos

    SQep = ((data[resp]-y_bar)**2).sum() #soma quadrática devido ao erro puro
    GLep = n-m #graus de liberdade devido ao erro puro
    MQep = SQep/GLep #média quadrática devido ao erro puro

    SQfaj = SQr - SQep #soma quadrática devido a falta de ajuste
    GLfaj = m-p #graus de liberdade devido a falta de ajuste
    MQfaj = SQfaj/GLfaj #média quadrática devido a falta de ajuste

    SQT = SQR+SQr #soma quadrática total
    GLT = n-1 #graus de liberdade total

    pve = SQR/SQT #porcentagem de variação explicada

    pmve = (SQT-SQep)/SQT #porcentagem de variação máxima explicável

    #imprime o modelo e a ANOVA
    print('Modelo:\n')
    for value,error,parametro in zip(coef,coef_var,entradas):
        if parametro == 'Intercept':
            print('{} = '.format(nome)+'{:.4e}±{:.5e} +'.format(value.item(0),1.96*error.item(0)**(1/2)))
        else:
            print('({:.4e}±{:.5e}) {} +'.format(value.item(0),1.96*error.item(0)**(1/2),parametro))
                  
    print("""
    Fonte de variação     Soma Quadrática      No de G. L.     Média Quadrática
    Regressão             {:>7.6e}              {:>2}              {:<6.6e}
    Resíduos              {:>7.6e}              {:>2}              {:<6.6e}
    Falta de ajuste       {:>7.6e}              {:>2}              {:<6.6e}
    Erro puro             {:>7.6e}              {:>2}              {:<6.6e}
    Total                 {:>7.6e}              {:>2}               
    % de variação explicada: {:.4f}
    % máxima de variação explicável: {:.4f}
    Teste F para ajuste do modelo: {:.3f} > {:.3f}
    """.format(SQR,GLR,MQR,
               SQr,GLr,MQr,
               SQfaj,GLfaj,MQfaj,
               SQep,GLep,MQep,
               SQT,GLT,
               pve,
               pmve,
               f.ppf(0.95,GLfaj,GLep),
               MQfaj/MQep,
              ))
In [26]:
fig = px.scatter(df, x='COF Médio', y="Corrente Média", trendline="ols")
fig.show()
In [27]:
results = px.get_trendline_results(fig)
print(results)
results.px_fit_results.iloc[0].summary()

                                      px_fit_results
0  <statsmodels.regression.linear_model.Regressio...

C:\Users\Caio\anaconda3\lib\site-packages\scipy\stats\stats.py:1604: UserWarning:

kurtosistest only valid for n>=20 ... continuing anyway, n=16
Out[27]:
OLS Regression Results
Dep. Variable: y R-squared: 0.358
Model: OLS Adj. R-squared: 0.312
Method: Least Squares F-statistic: 7.799
Date: Sat, 02 Jan 2021 Prob (F-statistic): 0.0144
Time: 15:56:56 Log-Likelihood: 161.97
No. Observations: 16 AIC: -319.9
Df Residuals: 14 BIC: -318.4
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const -6.061e-05 3.66e-05 -1.658 0.120 -0.000 1.78e-05
x1 0.0007 0.000 2.793 0.014 0.000 0.001
Omnibus: 0.529 Durbin-Watson: 2.615
Prob(Omnibus): 0.768 Jarque-Bera (JB): 0.483
Skew: -0.353 Prob(JB): 0.785
Kurtosis: 2.523 Cond. No. 92.6


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
In [28]:
anova(df,parametros,resposta,b,Vb,X,'Corrente')

Modelo:

Corrente = -6.0609e-05±6.33747e-05 +
(6.5522e-04±4.06665e-04) COF Médio +

    Fonte de variação     Soma Quadrática      No de G. L.     Média Quadrática
    Regressão             8.408273e-10               1              8.408273e-10
    Resíduos              1.509295e-09              14              1.078068e-10
    Falta de ajuste       1.509295e-09              14              1.078068e-10
    Erro puro             0.000000e+00               0              nan   
    Total                 2.350122e-09              15               
    % de variação explicada: 0.3578
    % máxima de variação explicável: 1.0000
    Teste F para ajuste do modelo: nan > nan
    

C:\Users\Caio\anaconda3\lib\site-packages\ipykernel_launcher.py:38: RuntimeWarning:

invalid value encountered in double_scalars
In [29]:
df['Prediction'] = b.item(0)+b.item(1)*df['COF Médio']
In [30]:
df['Residuals'] = df['Prediction']-df['Corrente Média']
In [31]:
fig = make_subplots(rows=2, cols=1,
                    row_heights=[0.25, 0.75],
                    shared_xaxes=True,
                    vertical_spacing=0.1,
                    subplot_titles=("Resíduos", "Dados",))
fig.add_trace(go.Scatter(mode = 'markers', x=df['COF Médio'], y=df['Residuals'], name = 'Resíduos'), row=1,col=1)
fig.add_trace(go.Scatter(mode = 'markers', x=df['COF Médio'], y=df['Corrente Média'], name = 'Pontos experimentais'),row=2,col=1)
fig.add_trace(go.Scatter(mode = 'lines', x=df['COF Médio'], y=df['Prediction'], name = 'Regressão'),row=2,col=1)

fig.update_xaxes(showgrid=False, row=1, col=1)
fig.update_yaxes(showgrid=False, row=1, col=1,tickformat='.0e')
fig.update_xaxes(title_text = 'Coeficiente de atrito (-)', row=2,col=1)
fig.update_yaxes(title_text = 'Corrente elétrica média (µA)', row=2,col=1,tickformat='.0e')

fig.update_layout()

fig.update_layout(
    title='Corrente elétrica vs Coeficiente de atrito cinético',
    showlegend=False,
    width=600,
    height=600
)
fig.show()